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Abstract 

We present a new approach to the problem of stationary viscoplastic 
duct flow as modelled by the Herschel-Bulkley model, with Bingham fluids 
included as a special case. 

While the mathematical formulation of this problem is conventionally 
based on a variational inequality, or equivalently, on a nonsmooth minimi¬ 
sation problem for the flow velocity, we suggest an alternative approach. 
Considering the Lagrangian dual in terms of the stress, rather than the 
velocity, turns out to be advantageous in numerous ways. The objec¬ 
tive functional possesses higher regularity, which ensures applicability of 
second order methods. Our numerical experiments with a trust-region 
SQP algorithm also demonstrate clearly superior performance compared 
to the widely used augmented Lagrangian method, although no artificial 
regularisation is introduced into the problem. 

Hence, besides providing a new theoretical angle to a classical prob¬ 
lem, our results also pave the way for an entirely new class of numerical 
approaches to simulating flows of viscoplastic fluids. 


1 Introduction 

Viscoplastic fluids are characterised by the existence of a yield stress, which 
marks the transition between viscous and plastic material behaviour. Examples 
of such materials are encountered in the consumer goods industry (toothpaste, 
paint), particularly in food processing (dough, tomato sauce) [I]. Viscoplastic 
models are also employed in geology and geophysics, for instance to describe 
flows of lava [5], lahars [3] or liquefied soil after a major earthquake [3]. 

The mathematical models named after Bingham [5] or Herschel and Bulkley 
[S] lead to problems of free boundary type, as a description of the interfaces 
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between liquid and solid phases is only implicitly contained in the governing 
equations. The nonsmooth transition between these two regimes and the non¬ 
uniqueness of the stress in plastic regions have always posed major obstacles for 
numerical simulations of viscoplastic fluids. 

Approaches to computational solutions of viscoplastic flow problems fall into 
two distinct categories. 

Methods that approximate the viscoplastic models with a purely viscous reg- 
ularisation were historically first developed. In this context we mention the work 
of Bercovier and Engelman [7] , the bi-viscosity model of Tanner and Milthorpe 
[5] and the widely used Papanastasiou regularisation [3]. More recently, de 
los Reyes and Gonzalez Andrade [13 E] employed Fenchel’s theory of duality 
and a Tikhonov regularisation of the dual problem to derive an algorithm that 
estimates the locations of yielded and unyielded regions as a by-product. 

All of these approaches share the advantage of fast convergence, since meth¬ 
ods of Newton type are directly applicable. They however become increasingly 
ill-conditioned the smaller the regularisation is chosen, i.e. the better viscoplas¬ 
ticity is approximated. Consequently, one must seek a compromise between 
stability and exactness of the method. Furthermore, replacing plasticity with 
a large, but not infinite viscosity means that methods of this first category are 
bound to fail with predicting circumstances under which a strictly viscoplastic 
flow would stop. This inadequacy with reflecting plasticity was demonstrated 
by Moyers-Gonzalez and Frigaard in m- 

Numerical methods from a second class solve the governing equations with¬ 
out prior regularisation. Duvaut and Lions HSllll] established a variational form 
of the Bingham flow problem in terms of a nonsmooth minimisation problem, 
or equivalently, a variational inequality. This pioneering work forms the basis of 
an augmented Lagrangian formulation of the problem and a corresponding nu¬ 
merical strategy commonly named ALG2. Originally, Glowinski m proposed 
this alternating direction method, which is widely used nowadays. For other 
techniques that avoid regularisation, such as further algorithms of Uzawa type 
or pseudo time relaxation, we refer to the review of Dean et al. m and the ref¬ 
erences therein. More recently, Aposporidis et al. m applied Picard iterations 
to a mixed formulation of the viscoplastic flow problem, which may or may not 
include regularisation. A few more unregularised algorithms are known [T81[T9] 
which, however, are only valid for a very limited range of applications. 

Not introducing any artificial regularisation to the problem means that plas¬ 
tic behaviour is exactly represented in the numerical solution, while the com¬ 
putation remains robust and well-conditioned. However, convergence tends to 
be significantly slower. In particular, methods like ALG2 and classical Uzawa 
algorithms do not achieve a quadratic or even superlinear rate of convergence. 
In addition, for methods that apply augmented Lagrangian techniques, the op¬ 
timal choice of free parameters in the algorithm, which strongly affect the speed 
of convergence, is still an open problem [^ pp I26-I27]. 

In our work, we follow a different approach. We first present a dual formu¬ 
lation for stationary Bingham or Herschel-Bulkley flow in terms of a linearly 
constrained minimisation problem. The objective functional possesses higher 
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regularity compared to the well-known classical functional of Duvaut and Li¬ 
ons. This feature opens up new possibilities for numerical solutions, by bring¬ 
ing second-order methods into play. In this paper, we choose to first discretise 
the optimisation problem with finite elements and propose a trust-region SQP 
method to tackle the resulting finite-dimensional problem. 

This paper is organised in the following manner: in Section [5] we introduce 
the model equations and a corresponding weak formulation. We also provide 
a summary of key results on the existence and uniqueness of solutions. Once 
we have presented the new, dual minimisation problem in Section [3] and its 
discretisation in Section 01 Section [S] is devoted to the numerical algorithm for 
its solution. Finally, we assess the performance of the new approach by carrying 
out a number of numerical experiments. 

2 Problem statement 

2.1 Strong formulation 

We consider the problem of flow through an infinitely long duct with homo¬ 
geneous cross-section. We assume the cross-sectional domain 17 C to be 
Lipschitz with boundary 917 = T. Furthermore, we let the boundary be decom¬ 
posed into measurable, disjoint Dirichlet and Neumann sections F = Fd U Fn, 
where Fd is required to possess positive measure. 

The Bingham and Herschel-Bulkley models provide constitutive relations 
between the stress r and the rate of strain 7 ; both are vector fields 17 — 

For the case of duct flow considered here, the rate of strain is given by Vy, the 
gradient of the scalar flow velocity y : 17 ^ K. through the cross-section. Model 
parameters include a plasticity threshold or yield stress tq > 0, for Bingham 
fluids additionally a plastic viscosity parameter /r > 0 and for Herschel-Bulkley 
fluids a consistency k > 0 as well as a power-law exponent a > 1. 

The flow is driven by a pressure gradient or volumetric force density / : 17 —>■ 
R. 

In the classical strong formulation, stationary duct flow of a Herschel-Bulkley 
fluid is governed by the system 


1 • 1 a —2 • 1 "y 

t = k| 7| 7 + To-p^ 

if 7 ^ 0 

( 2 . 1 a) 

|t| < To 

if 7 = 0 

( 2 . 1 b) 

with conservation of momentum 



—div T = f 

in 17 

( 2 . 1 c) 

and either no slip or perfect slip (symmetry) boundary conditions 


y = 0 

on Fd 

( 2 . 1 d) 

T ■ fi = 0 

on Fn- 

( 2 . 1 e) 
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Here n denotes the outward normal unit vector on the boundary section Fn- 
The special case of a Bingham fluid is recovered through the choices a = 2 
and jj. = K. (I2.1al) and (I2.1bl) then assume the simplified form 

X = /r 7 + To ^ if 7 ^ 0 (2.2a) 

l7l 

|t| < To if 7 = 0. (2.2b) 

2.2 Weak formulation 

A weak formulation of the Bingham flow problem is introduced in [mE]. An 

extension to general Herschel-Bulkley fluids is given in m- Similar to these 

works, we consider the following function spaces: 

For Lebesgue spaces of order a we apply the standard notation L“(B). 

Spaces of R^-valued functions are set in boldface. Equipped with the canonical 

norms IHI^, T“(B) and L°‘{n) are Banach spaces. If a = 2, i.e. within the 

Bingham scenario, we obtain Hilbert space structure with respect to the or 

scalar product (•,•)■ The pairing between or L°‘{n) and their dual 

/ / 

spaces L°‘ (B) or L°‘ (B), respectively, is referred to as (•, •)^, where the dual 
index a' is defined through l/a + l/a' = l. As usual, dual spaces are marked 
with an asterisk. 

With IF^’“(B) we denote the Sobolev space of L“(B)-functions with first 
derivative in L“(B). 

For the Euclidean scalar product in R'^, we use the notation •, for the a-norm 
in we use | • |a- If a = 2, we may omit the subscript. The symbol C is used 
as a generic constant. 

Furthermore, we introduce the subspace of admissible velocity fields and 
velocity test functions 

y := {y e IFi’“(B) : ylr^, = O} C IF1’“(B) 

For the ease of notation, we define the bi-linear, if a = 2, otherwise non-linear 
form a : y X y ^ R 

a(y, z) := k Vy, 

and the non-linear functionals j :Y —>■ R 

j{y) ■= tq J |Vy| dx 

n 


as well as c : y —R 


c(y) := -a(y,y) = - f |Vy|“ dx. 
a a J 


We justify that with the above definitions, a and j are indeed well-defined. 
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Remark 2.1. (a) Since y G we have Vy G L°‘{Vt) and thus 


|Vyr-"Vy 


n 

< C 


iVyir'Vy 


iVj/ir'Vy 


da; 


da; 


= cj |Vy|“' da; 

n 


= CJ IV2/I2 

n 




da; 




Consequently, |V 2 /|“ ^ Vy G i“ (fl). Hence, the definition of a is sensible. 

(b) Since fl has finite measure, L“(H) C L^(r2). In particular, |Vy| G L^(n) 
which guarantees that the integral in j is well-posed. 

Definition 2.2. Let f & L°^ (H) be a given force density. We call a velocity 
field y €Y a weak solution to the Herschel-Bulkley problem if 


a{y,z-y)+j{z)-j{y)>{f,z-y)^, ^zeY. (2.3) 

It turns out that this variational inequality of the second kind is equivalent 
to the first order necessary optimality condition for the minimisation problem 


inil{y)=c{y)+j{y)-{f,y)^, (2.4) 

yeY 

cf. [21], [221 p 8]. 

Huilgol and You [21] conclude that “there are gaps as far as the existence and 
uniqueness of a solution to the pipe flows of the three fluids under consideration” 
[H p 141] (the authors consider Bingham, Herschel-Bulkley and Casson fluids). 
In fact, it turns out that both existence and uniqueness can be established with 
well-known tools from convex optimisation; 

Theorem 2.3. Problems (]2.3I) and (12.4p are equivalent and possess a unique 
solution y* GY . 

Proof. Existence. A standard proof for the existence of a solution to (12.41) can 
be found e.g. in Theorem 1.1 and Remark 1.2 of [23l pp 7-8]. It requires weak 
lower semi-continuity of the objective / and 

I{y)^+oo as ||?/||wi>-(n)-)> 00 (2.5) 

to ensure boundedness of any minimising sequence. 
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The first assumption follows immediately since I is both continuous and 
convex, cf. Theorem 3.3.3 in p 93]. 

To verify the second assumption, we use the following generalisation of 
Poincare’s inequality: there exists a constant C > 0 independent oi y € Y 
such that 





( 2 . 6 ) 


A proof for the special case a = 2 is given in p 355]. An analysis shows 
that it remains valid for 1 < a < oo. 

Since y = 0 almost everywhere on Td, the last term in (12.61) vanishes. Con¬ 
sequently, using (12. 6 |) along with Holder’s inequality, we obtain 



1 V?/|“ dx + To J 


|Vyl dx - {f,y)a 


^ -I- 0 — ll/IU' ||y||a 

^ C (ll2/llivi><»(n) — ||2/||wua(n)^ 
—>■ -boo as |] 2 /|lvvu<»(n) ^ oo. 


Uniqueness. As a consequence of the strict convexity of /, a solution of (12.4|) 
is unique. 

Equivalence. We refer to Theorem 1.6 and Remark 1.12 in [23l pp 12-13]. □ 


2.3 Constrained formulation 

Augmented Lagrangian methods are a common choice for the numerical solution 
of nonsmooth optimisation problems. One first introduces a new variable q for 
the rate of strain Vy, where q and Vy are linked in a constraint. This trick 
allows to decouple nonsmoothness and nonlinearity on the one hand from the 
linear velocity term on the other hand. The constraint is then relaxed and 
penalised in an augmented Lagrangian. 

The unconstrained problem (|2.4I) rewritten as artificially constrained prob¬ 
lem becomes 

, min = - /kr d 2 :-bro /"iq] dx-(/,y)a (2.7a) 

(y,q)eyxi“(n) aj J 

n n 

subject to q — Vy = 0. (2.7b) 

We define the corresponding Lagrangian L :Y x x L°‘ (H) —)■ K by 

L{y, q,T) = J \q\°‘ dx -b To J \q\ dx - (/, y)a - (r, q - Vy)^. (2.8) 

n n 

We will now answer the question for the existence of Lagrange multipliers. 


6 


Proposition 2.4. Let {y*,q*) £ Y x L°‘{Q) be the solution of (12.7|) . There 
exists a Lagrange multiplier r* G Z/“ (LI) such that the KKT conditions 


DyL{y*,q*,T*)=0 (2.9a) 

dgLiy*,q*,T*)3 0 (2.9b) 

D^L{y\q\T*)=0 (2.9c) 

hold, i.e. 

(r*,Vz)„ = (/,z)„ V^GF 

(2.10a) 

— r*, r — q*)^ + To J |r|da; —tq J |q*| da; > 0 Wr £ L°‘{Ll) 

(2.10b) 

q* = Vy* a.e. in Ll. 

(2.10c) 


Here, Dy denotes the Frechet derivative with respect to y and dq the subdijfer- 
ential with respect to q. 


Proof. The assertion is an immediate consequence of the fact that Slater’s con¬ 
straint qualification (SCQ) is trivially satisfied here (cf [551 Thm. 3.34]). □ 

Physically, a Lagrange multiplier r associated with the constraint (I2.7b|) 
can be interpreted as an admissible stress m- Therefore, the (Lagrangian) 
dual problem to (I2J1) can be seen as an optimisation problem in terms of the 
stress T instead of the veloctiy y. We will now follow this dual approach. 


3 Dual formulation 

3.1 The Lagrangian dual 

Lagrangian duality expresses the fact that a constrained optimisation problem 
can be represented in two distinct, but often equivalent ways: the primal prob¬ 
lem in terms of the primal variables or the dual problem in terms of the dual 
variables (Lagrange multipliers). 

From the Lagrangian L in (12.811 . one obtains a primal function Lp by max¬ 
imising in the dual variable 

Lp{y,q)-= sup L{y,q,T) 

-reL°‘'{n) 

and a dual function L^ by minimising in the primal variables 


Ldir) 


inf 

{y,<i)eY xi“(n) 


Liy,q,r). 
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The primal and dual functions are permitted to assume values on the extended 
real line, K U {+00} and K U {—00}, respectively. Then the primal problem is 
given by 


min 

(y,q)GYxL‘-(n) 


Lp{y,q) 


and the dual problem by 


max Ld('r). 


Since we have just seen that SCQ holds, we actually have strong duality 


min 

(y,q)GVxZ,“(n) 


Lp{y,q) 


max Ld(T), 


i.e. there is no duality gap. Furthermore, with solutions {y*,q*) and r* of 
the primal and dual problem, respectively, {y*,q*,T*) is a saddle point of the 
Lagrangian (primal-dual solution), see e.g. [IHl Thms. 4.7, 4.9]: 

L{y\q*,T) < L{y\q*,r*) < L{y,q,r*) y{y,q,T) eYx x L^'{n). 

Calculating the primal function from the Lagrangian L in (12.8|) yields 


Lp{y,q) 


a J 1^1 

n 


- 1-00 


if q^Vy 
otherwise 


Clearly, the primal problem minLp(?/, q) is equivalent to (12.71) . 

/ 

Let us now turn to the dual problem. We assume t G L°' (n) is given and 
we will now calculate the corresponding value of the dual function Ld(T)- 

Initially we observe, since L depends linearly on y, that L^{t) = —00 if 
(r, V-)a (/, ■)a in Y*. This is a weak formulation of —divr ^ /, stating that 
T violates conservation of momentum (I2.1cp . 

Otherwise, if (t,V-)q, = {fr)a holds as equation in Y*, we have to find 
q* G L°‘{n) which minimises the remaining convex functional 


q = 


argmin- / \q\°‘dx + tq / |q| da; - (t, q)„. (3.1) 

qeL‘‘(n) ct J J 


(r2) 

Q, Q, 

Lemma 3.1. The unique solution of dSH) is 


q = 


1 


.l/(a-l) 


Nl/(a-l) "T 
-To)+ -j—r. 


(3.2) 


where (•)+ := max(0, •) represents the hard-thresholding operator for truncation 
below zero. 





Proof. Feasibility. First of all, we show that the above q* is admissible, i.e. 
belongs to Since t £ L°‘ (O), we have |t| — tq G and thus 


kl -T-o)+“ 


a 

< C 


( 11 - 12 - 7 - 0 )+“ 


F |2 


dec 


|x|2-ro)r ^ 


1-2 


dx 


= C/(|x|2-ro)r ^dx 


= C / (|x| -ro)+ dx 


< C / ||x|—To|“ dx < 00 . 


Consequently, (|t| - to)+ ^ G L°^{Vl). 

Optimality. In analogy to Proposition 12.41 we conclude that the minimising 
argument satisfies the first-order optimality condition 


(K\q*r^q* 


— T,r 


— q*)^+To j |r| dx—To y |q*| dx > 0 Vr G L°‘{^1). (3.3) 


Due to the strict convexity of the problem, this condition provides a character¬ 
isation for the unique solution q*. 

We set s := r — q*. Using the above q* we obtain for the first term 


iK\q 


= ( (|t| -To)+|^ -T,s^ 

J (l'^|-'^o)|^-sdx-(T,s), 


IItUto} 


/ T • S dx — To / ■;—: ■ S dx. 

J J fI 

{|-rl<To} {|-r|>i-o} 


These two integrals can be estimated from below: from the Cauchy-Schwarz 
inequality we derive 


T ■ s dx < 


|t||s| dx < To 


\s\ dx = To 


{Ul<T-o} 


{UI<T-o} 


{|T|<ro} 


J |s-hq*|dx. 


{UI<T-o} 


(3.4) 
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since q* = 0 almost everywhere on {ItI < tq}. Similarly we conclude 


To 


: • S dx = To 


■j—7 ■ s dx + To 

X 


|q*| dx - To / |q*| dx 


{|t|>to} 


{|t|>to} 


{kl>To} 


= To 


{|t|>to} 


< To 


{|-r|>To} 


= To 


{k^T-ol 


dx-roy |q*|dx 
n 

^ + dx-To J |q*|dx 

n 

s + q*| dx-To J |q*|dx, 


where we used q*/\q*\ = t/|t| almost everywhere on {|t| > tq}, which proves 
the assertion. □ 

By substituting this solution q* in the Lagrangian (12.81) and simplifying, we 
obtain the dual function 

, /(|t| - To)X dx if (t, V-)a = (/, ■)a in Y* 

Lo{t) = { n 

oo otherwise 

Instead of the dual problem max Ld(T), we will from now on consider the equiv¬ 
alent problem 


At) = J (1^1 - ^o)+ 

a 

subject to (t, V-)a = (/, •)« in Y* 


min 

xei“'(n) 


(3.5a) 

(3.5b) 


and, with a slight abuse of terminology, refer to p.5l) as the dual problem. 

The numerical method we propose will be based on this characterisation of 

Herschel-Bulkey duct flow. We remark that (13.5p is convex, but generally not 

/ 

strictly convex. Hence, while a solution r* G L°‘ (H) is always guaranteed to 
exist, it needs not be unique. 


3.2 Optimality conditions 

We already saw, since SCQ holds for the primal problem, that the primal and 
dual problem possess solutions and that strong duality holds, i.e. both problems 
are equivalent. We could also arrive at this conclusion by verifying, say, the 
linear independence constraint qualification (LICQ) for the dual problem, which 
for any dual solution r* S (fl) implies the existence of a (unique) Lagrange 
multiplier y* G Y. 


10 






Due to the convexity of problem (13.51) , the corresponding KKT conditions 


1 


K 




(kl -'ro)+ ' j -((T,V?/)„ = 0 Vcr e (D) (3.6a) 

= (/> 2)a e Y. (3.6b) 


are necessary and sufficient for optimality. This justifies the following definition; 

Definition 3.2. We call a pair (x, y) G (D) x T a weak solution to the mixed 
Herschel-Bulkley problem, if they satisfy (ED). 

One can relate (13.61) back to the strong formulation (12.11) we started with; 
looking at (I3.6ap . a corresponding strong form reads 


In plain words, the rate of strain can be recovered from the stress by a soft- 
thresholding operation, combined with a power law if a ^ 2. By rearranging 
for the stress r, we see that this is equivalent to the viscoplastic constitutive 
relations of the Herschel-Bulkley model in (I2.1al) - (l2.1bl) . cf p71128) ; 


Vy = 


1 

Kl/(a-l) 


(|t| - 'To) 




T = K|Vy|“ ^ Vy-bToi||| 

kl < tq 


iiVy^O 
if Vy = 0. 


We already identified (I3.6bl) as a weak formulation of conservation of momen¬ 
tum (j2.1c|) . The boundary conditions (j2.1dp - p.le|) are incorporated through an 
appropriate choice of function spaces for the solution and the test functions. 


3.3 Regularity of the objective 

Numerical optimisation algorithms typically benefit in terms of their speed of 
convergence when higher derivatives of the objective and any constraints exist 
and are used. Due to the nonsmoothness in j, the conventional functional I in 
(1^ is not Frechet differentiable whenever a vanishing velocity gradient occurs. 
In contrast, the stress functional J in (13.51) is generally continuously Frechet 
differentiable in t. Depending on the value of a, additional regularity may be 
achieved. 

For shear-thickening Herschel-Bulkley fluids, a > 2, no higher Frechet deriva¬ 
tives of the objective J are available. This is due to the fact that the optimality 
system (13.61) is not differentiable at the interface between viscous and plastic 
regions; 

Let T* G T“(H) with |t»| = tq be the limit of a sequence (T")„gN with 
|r"| < To for all n G N and |r”| ^ tq as n ^ oo. Let (cr”)„gf^ be another 
sequence with lim„_>oo o'n = t*, \cr'^\ > tq for all n G N and |cr”| \ tq as 
n —>■ 00. J is twice Frechet differentiable at every r" and every cr”. From 
(|3.5ap we observe, however, that 

hm ||J"(r")||=0 


II 





















while 


lim ||J"(c 


= oo. 


This sudden jump from an infinite to a vanishing viscosity is not a flaw of the 
analytical framework chosen here, but rather an unphysical feature of the model. 

For shear-thinning Herschel-Bulkley fluids, 1 < a < 2, one can conclude 
from (|3.5a|) that J is at least twice continuously Frechet differentiable. 

For Bingham fluids, a = 2, J is once, but not twice continuously Frechet dif¬ 
ferentiable. However, J' still satisfies a slightly weaker notion of differentiability 
that is often applied in the context of semismooth Newton methods. 

Definition 3.3. Let K^L be Banach spaces, D C K an open subset. 

A function F : K D D ^ L is called Newton differentiable (or slantly 
differentiable) in the open subset U C D, if there exists a family of mappings 
G :U ^ C{K,L) such that 


lim 


||F(a; + h)- F{x) - G{x + h)h\\^ 


= 0 , 


Vx e U. 


(3.7) 


For a proof that J' is indeed Newton differentiable, we refer to Lemma 3.1 

in P5] . 

This result has far-reaching consequences: Sequential Quadratic Program¬ 
ming (SQP) methods attempt to solve the system of first order necessary opti¬ 
mality conditions with methods of Newton type. Since, independent of a, the 
functional / does not possess a continuous Frechet derivative where Vy = 0, 
in particular no second Frechet or Newton derivative, SQP methods are not 
applicable to this problem. With the alternative approach suggested here that 
minimises J, Newton differentiability (or even higher regularity) of the optimal¬ 
ity system is warranted for the most important applications with shear-thinning 
and Bingham fluids, 1 < a < 2. Therefore, the new approach to formulating 
viscoplastic flow allows us to use powerful numerical algorithms from the SQP 
framework. 

A similar duality approach was chosen by de los Reyes and Gonzalez Andrade 
nnmi] for Bingham fluids. However, the authors used a different splitting in 
the primal objective, which results in a more complicated dual problem. They 
then introduce a Tikhonov regularisation into the dual objective and use a 
semismooth Newton method to solve the resulting optimality system. We could 
certainly follow an analogous procedure here and solve a penalised dual problem 
with a semismooth Newton method for Bingham flow, and a smooth Newton 
method for shear-thinning Herschel-Bulkley flow. However, the very simple 
structure of the dual (13.511 allows us to even solve the unregularised problem. In 
our numerical method, this will require the application of trust regions, resulting 
in exactly viscoplastic solutions. 
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4 Discretisation 

For the numerical solution of (1^ . we follow a first-discretise-then-optimise 
approach. That is, we first discretise both the objective functional and the 
constraint in order to obtain an optimisation problem that is posed in finite 
dimensions. We are then in the position to apply well-established methods to 
solve this problem numerically. 

Let be a polygonal domain that approximates 17, 0 7^ F^ C the set 
of all edge segments that correspond to Fd, analogously F^ = 917^ \ Fq and 
a regular triangulation of the closure 17^. We discretise with Pl-PO finite 
elements, i.e. discretise the stress with piecewise constant functions, the velocity 
with piecewise linear functions. Denoting with the space of polynomials in 
with degree k or less, we therefore introduce the following finite-dimensional 
spaces: 

:= jr'* = e X“'(17'*) : t^\t £ Pq, VT £ T^,j £ {1,2}} C i“'(17'*) 

:= £ (7(17^) : /|t £ Pi, VP £ T'*} C ^^’“(17''). 

Furthermore, we define the subspace of discrete velocity functions 
Y’^ :={y'‘£P'‘:/ = 0on F^} C 
A basis of the space is given by 


X* 


XfeCi if i = 2/c - 1 

Xke2 iii = 2k 


1,..., 2nT- 


(4.1) 


Here, Xk is the characteristic function of the triangle Tk £ is the number 

of triangles in and e}, 6*2 are the canonical unit vectors of . 

In the following we will assume that the un triangle nodes {x^ are sorted 
such that the first ni nodes, nj < un, lie in the interior or on the Neumann 
boundary of the domain, {x^ )^Li C 17^ U F(^, the remaining ud = un — nj nodes 
on the Dirichlet boundary, c F^. By (()j, j = 1, ...,nN we denote the 

hat function of the node xP Then form a basis for and {(t>j)^Li span 

the subspace Y^. 

Using these spaces, a discrete counterpart of (13.51) reads as follows: find a 
solution for the minimisation problem 


n 

subject to (r^ V-)^ = (/, •)„ in {Y’^)* (4.2b) 

An equivalent representation of (14.2|) can be obtained through a standard 
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(4.3a) 


procedure, which yields an optimisation problem posed in 

fc=l 

subject to D^t = (4.3b) 

Dh g ]graix2nT^ g ^ g ]^2raT gj-g ^]gg representations of (•, V-)q,, 

{fr)a or in terms of {Xi)'i2i or {4>j)'jLn respectively. We use fk as an 
abbreviation for (T2 fc-i,G K^. 

Proposition 4.1. Problem (14.31) has a solution f* G 

Proof. First, the admissible set for potential minimisers r G as defined 

through the constraint D^t = is nonempty: since, due to the inf-sup stable 
definition of the finite element spaces, has full rank ni [21] and is 

invertible. One particular feasible point is given by 

T := {D'^Y{D^{D’^Y)-^f^. 

Furthermore, is convex, continuous and -l-oo as |t| —>• oo. 

In analogy to Theorem 12.31 we conclude that there exists a minimiser r G 

□ 


Proposition 4.2. If r* G is a solution to (14.3F then there exists a 

Lagrange multiplier y* G such that (r*, j/*) satisfies the KKT conditions 

Vj'^iT*) - {D'^)^y* =0 (4.4a) 

jjhf* ^ (-4 4 ^) 


where VJ^{t*) = given by 


V,G'‘(fl = ;^(|TV|-ro)r ^ 


(4.5) 


Proof. This result follows immediately from the inf-sup stability of the chosen 
finite element spaces, i.e. is onto. Hence, LICQ holds for any r G R"'^. □ 


5 Trust-region SQP algorithm 


We will now present the trust-region SQP algorithm that we propose to solve 
problem (14.3F It is based on the Byrd-Omojokun method [301 EH using an 
implementation similar to [32] . 

In order to guarantee sufficient smoothness of the objective, we assume from 
now on I < a < 2. Then the Hessian J^(t) is locally bounded. From (14.51) 
we obtain the following explicit representation: 




'ir) = 


\Tk\ {\Tk\-To)l 

^(a-I) 


-H{Tk), 


(5.1) 


14 





where the symmetric 2x2 matrix H has the components 

/ill = (a - {\fk\- To)+ + Tffe_i|ffc| 

hi2 = -T2k-lT2k ((a - 1) (|Tfe| - ro)+ - |Tfc|) 

h22 = (a - l)/ffc-l (|Tfe| - To)+ + 7ffc|Tfe|. 

5.1 Sequential Quadratic Programming 

The rationale behind Sequential Quadratic Programming (SQP) is to approxi¬ 
mate the nonlinear minimisation problem (14.31) with a sequence of linear-quadratic 
problems, i.e. optimisation problems with a quadratic objective and linear con¬ 
straints. For problems with equality constraints only, like here, a basic SQP 
method is equivalent to applying Newton’s method to the first order optimality 
conditions (lia . Given an iterate ,y^) ik = 0,1,2,...), Newton’s method 
attempts to find (dr, 5y) such that 

- D^Sy = (5.2a) 
DSt =- f). (5.2b) 

To simplify the notation, we have suppressed the superscript h and we will 
continue to do so from now on. A superscript k denotes a function evaluation 

at T . 

The system (lOl) can be identified with the KKT conditions of a linear- 
quadratic problem. Indeed, one could equivalently try to obtain a step dr by 
solving 

mindx VJ^-|--dx V^J^dx (5.3a) 

subject to D6t -(- Df^ — / = 0. (5.3b) 

The linearity of the constraint (I4.3bl) allows us to implement some simplifi¬ 
cations to the general numerical method. Given any x € we can find a 

projection x° onto the manifold {x S : Dr = /} by setting 

x° := X - A-^D^{DA-^D^)-^{Df - /), (5.4) 

where A is a diagonal matrix that contains the triangle areas 

A := diagdTil, |Ti|,..., |T„.j,|, |r„.j,|). 

The matrix DA~^D^ is the usual discretisation of the Laplacian with piecewise 
linear Lagrange elements and A~^D^ is a discrete analogue of the gradient. 
With (15.4|) . the constraint (I5.3bl) simplifies to 

DSt = 0 (5.5) 

and all iterates x^, k = 0,1,2,... are automatically feasible with respect to 
gjb]). 
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Let G R"'^ x M”' be an exact solution of the KKT conditions (14.41) . 

Given an iterate G R"'^, the overdetermined system (I4.4all may not have 
a corresponding solution G R"L However, we can obtain a solution in the 
least-square sense by solving 

f ■- argmin \Wj’^ - D^y\^ = {DD^)-^DVj'^. (5.6) 

Later on, we will use (15.61) to get estimates of y* from the iterates r^. 

5.2 Trust region 

Trust-region methods impose an additional constraint on an optimisation prob¬ 
lem by limiting the step size with a trust radius > 0. We apply the 

trust-region concept to the subproblem (15..'ll) with (15.51) . Intuitively, one could 
say that we only trust the Taylor approximation of the objective in (15.31) to be a 
sufficiently accurate representation of the exact objective within a ball of radius 
around r^. 


min^T VJ^ -I- -<5 t 

(5.7a) 

subject to D5t = 0 

(5.7b) 

< A^r. 

(5.7c) 


Through this amendment, we enforce a finite step length even if the Hessian is 
not positive definite. From (EU, we observe that zero eigenvalues will occur 
whenever \Ti\ < tq on a triangle T, (i G {1,..., ut})- Hence, we can only assume 
the Hessian J to be positive semidefinite, but not positive definite. 

Generally, trust-region SQP methods require the solution of two subprob¬ 
lems, often referred to as the vertical (or normal) and the horizontal (or tangen¬ 
tial) subproblems. With an affine constraint like (I5.3bp . the set of admissible 5t 
that satisfy both this constraint and the trust-region constraint may be empty. 
Then, the purpose of the auxiliary vertical subproblem is to replace the equal¬ 
ity constraint with a suitable relaxation. However, due to the initial projection 
that we suggest, (I5.7bl) and (I5.7cl) are clearly not mutually exclusive (5 t = 0, 
for instance, satisfies both constraints). As a result, for the problem under con¬ 
sideration, we may do without the vertical subproblem. Without any further 
modifications, the horizontal subproblem is given by (E3. 

5.3 Solution of the horizontal subproblem 

We are looking for an efficient method that provides an approximate solution to 
the horizontal subproblem (15.71) . While the classical Gonjugate Gradient (GG) 
method assumes a positive definite Hessian, Steihaug [33] and Toint [33] devel¬ 
oped an extension that remains valid for positive semidefinite or even indefinite 
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matrices. Our implementation of the projected CG-Steihaug algorithm is an 
adaptation of the algorithms in [32l [35] . 

Algorithm: Projected CG-Steihaug (CGS) 

Data: k € No, divtol > 0, abstol > 0, reltol > 0, > 0 

Set z° = 0, = VJ^ / = - D^{DD^)-^Dr^, d = 

if < abstol then return = 0 for j = 0, 1 , 2 , ... do 

if (f < divtol then 

/* Small or zero curvature */ 

/* Intersect step with trust region boundary */ 

Find s > 0 such that \z^ + sd’ p = A^j^; 
while Armijo condition (15.8|) violated do 
I Backtracking s —>■ s/2 

end 

return + sd' ; 

end 

Compute z^ = z^ + / (S' ; 

if |£^| > A^pj^ then 

/* Step exceeds trust region */ 

/* Intersect step with trust region boundary */ 

Find s > 0 such that \z^ + scS p = A^j^; 

return = z^ -|- s(S ; 
end 

Set 5^+^ = -b 

^i+l = pJ + i _ £)Tp^T)-l£,-J + l. 

if \J < reltoly^9°^?° then 

/* Solution has converged */ 

return St^ = 
end 

Set / (g^^r^'^] 

Compute = —g^~^'' + fi^'^^d! ; 

end 

The Armijo condition 

Q{z^ + sd') < Q(z^) + 7sVg(z^yd^' (5.8) 

with 7 = le-2 guarantees sufficient decrease in the model objective Q{z) := 
z VJ^ -b iz 'S/^J^z even if very small curvature is encountered. 

One can show that while the model objective (I5.7al) decreases monotonically 
with every iteration, the iterates strictly increase in their norms, > 15-^1 
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[551 p. 172]. These properties justify the stopping criteria of the first two if 
statements within the for loop. 

5.4 Update of the trust radius 

To analyse the quality of a step St computed by Algorithm CGS, we compare 
the actual reduction in the value of the objective 

ared := j’^ - j'^+^ (5.9) 

with the reduction predicted by the quadratic model (I5.7al) 

pred := -St V- -St V^J^St . (5.10) 

If the actual reduction is sufficiently large, ared > rj pred with some parameter 
r] e (0,1), then the step is accepted. A particularly good agreement of ared 
with pred indicates that we may even extend the trust region in the next step. 

Otherwise, if ared < r] pred, we infer that the discrepancy between the 
exact objective and its Taylor approximation is too large. Consequently, we 
reject the step computed in Algorithm CGS and shrink the trust region. 
Following [22], we update the trust radius as described below. 

Algorithm: Update of the trust radius (UTR) 

Data: A: e No, ry G (0,1), > 0, > 0 

Compute p = ; 

if p > T] then 

/* Step is accepted */ 

Set St = St ■, 

/* Increase trust radius depending on model accuracy */ 
if p > 0.9 then 

I = min(max(10|dx |, A^p,), A™“); 

else if p > 0.3 then 

I = min(max(2|dT^|, A^jp), A““); 

end 

else 


/* Step is rejected 

*/ 

Set St = 0; 


/* Decrease trust radius 

*/ 

A^"^^ = max(0.1, min(0.5, ^^)) dx |; 



end 
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5.5 Algorithm TRS 

All in all, the computation of a solution to (14.31) follows the iterative scheme as 
outlined below. 

Algorithm: Trust-region SQP (TRS) 

Data: abstol, reltol > 0,r] £ (0,1), > 0, G (0, A™“] 

Choose f G ; 

Compute projection t° with (15.41) : 
for fc = 0,1, 2,... do 

Evaluate VJ^ and compute with (15.61) : 

if |VJ^ — < abstol and A reltol|^| then 

I return r , y 
end 

Evaluate V^J^; 

Solve the tangential subproblem (15.7p for St with Algorithm CCS; 

Compute T~^ := -b 
Evaluate and J~^; 

Compute ared and pred from (15.9|) and (15.101) : 

Use Algorithm UTR to accept or reject 

end 

We note that a stopping criterion which measures the relative difference 
between subsequent iterates of the stress is not appropriate, since in general the 
stress is not uniquely determined. 


6 Numerical experiments 

We will now present numerical results to demonstrate the performance of a 
MATLAB implementation of Algorithm TRS. 

In all examples, we initialise the algorithm with r = 0, i.e. the initial iterate 
is given by r° = — A~^{DA~^)~^f according to ()5.4I) . Eurthermore, we 
set abstol = reltol = le-4, divtol = le-10, A!pj^ = 10, Aip“ = le5 and y = 
le-1. Eor the problem parameters we assume y = k = 1 and / = 1. 

Numerical solutions for Bingham and Herschel-Bulkley flows through various 
duct geometries are available in the literature umi]. These were computed 
with Algorithm ALG2, an Uzawa type method based on an augmented La- 
grangian formulation of the nonsmooth minimisation problem (12.41) . As [21) 
already contains a discussion of the results in [12] , we will use the data provided 
in [21] as reference solutions. 

In order to compare the runtimes of Algorithms TRS and ALG2, we also 
provide results which we obtained with our own MATLAB implementation of Al¬ 
gorithm ALG2. Here we use the discrete analogue of the augmented Lagrangian 

Lriy, q, t) = c{y) + j{y) - {/, y)a -{T,q- + ^ J l^y- qp da;, (6.1) 

n 
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discretised with the same finite elements as for Algorithm TRS. We remark 
that for a < 2 it is not obvious whether the penalty term in (lO) is well-defined 
since Vy and q do not necessarily belong to This problem dissolves in 

the discrete setting. Alternatively one might consider a non-quadratic penalty 
term, e.g. by replacing the exponent 2 with a. In that case, however, the 
penalty term will not give rise to a Laplacian in the optimality conditions, but 
rather a quasilinear elliptic operator similar to an a-Laplacian. Consequently, 
ALG2 would require the solution of an additional nonlinear equation in every 
single iteration. 

We set r = 10, which by trial and error appears to be an optimal choice 
to minimise the number of iterations. To ensure consistency with the set up 
of Algorithm TRS, we terminate Algorithm ALG2 as soon as the residual of 
the first-order optimality conditions measured in the oo-norm no longer exceeds 
abstol. We also require subsequent iterates of the discretised velocity y and 
the discretised rate of strain q to have a relative difference of at most reltol. 
For Newton’s method, which is required to find a new iterate for q if a 7^ 2 (cf 
EH), we also use abstol and reltol as absolute and relative tolerances. 

Our meshes are generated with the M ATLAB built-in Delaunay triangulation. 
The programs are executed with MATLAB R2013a 64-bit on a laptop with an 
lntel®CoreT'^i5 CPU 4x2.53 GHz. 


6.1 Flow through a cylindrical pipe 

We consider the classical test problem of flow through a cylindrical pipe with 
radius 1 and homogeneous Dirichlet boundary conditions. A typical solution is 
depicted in Figure [H The graph illustrates very well the existence of a plastic 
region with stagnant flow which is separated from a region with regular viscous 
flow. 

For this quasi one-dimensional problem, an analytical solution is known. 
With R := |a;|, i?o := 2ro and /3 := l/(a — 1) it reads [36] 


y{R) = 


-(R-Roy+P) 


0<R<Ro 
Ro< R<1' 


( 6 . 2 ) 


Huilgol and You [21] provide results which they computed with their im¬ 
plementation of Algorithm ALG2 for this problem. They consider the cases of 
a = 2, a = 1.75 and a = 1.5 on three different grids and with tq = 0.1 or 
To = 0.2. The authors write n for a — 1, Od for 2ro and u/U for y. 

While Huilgol and You report that a “convergence tolerance e = 10“^ is kept 
fixed” jUJ p 130], it unfortunately remains unclear what this tolerance refers 
to. Their numerical results suggest that their stopping criterion only measures 
a relative difference between subsequent iterates of the flow velocity. In this 
case, their implementation of ALG2 would neglect to check how well the first 
order necessary optimality conditions of the minimisation problem are satisfied. 

In Tables |T] to |3| we present some typical results of Algorithm TRS and com¬ 
pare them with those of our implementation of Algorithm ALG2 as well as the 
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To 

UN 

Relative error 
TRS ALG2 

Iterations 
TRS ALG2 

CPU time (s) 
TRS ALG2 

Speedup 

0.1 

559 

1.48e-3 

2.27e-3 

3 

73 

0.06 

0.36 

6.0 

0.1 

1129 

6.35e-4 

1.42e-3 

2 

73 

0.06 

0.66 

11 

0.1 

2169 

3.40e-4 

1.24e-3 

1 

73 

0.07 

1.19 

17 

0.2 

559 

2.19e-3 

3.01e-3 

58 

73 

0.49 

0.40 

0.82 

0.2 

1129 

8.97e-4 

1.75e-3 

10 

73 

0.16 

0.66 

4.1 

0.2 

2169 

5.30e-4 

1.41e-3 

2 

73 

0.09 

1.18 

13 


Table 1: Bingham flow through the cylindrical pipe, a = 2, using un triangle 
nodes. 


To 

Un 

Relative error 
TRS ALG2 

Iterations 
TRS ALG2 

CPU time (s) 
TRS ALG2 

Speedup 

0.1 

559 

1.97e-3 

2.67e-3 

11 

67 

0.14 

0.76 

5.4 

0.1 

1129 

8.79e-4 

1.55e-3 

2 

67 

0.07 

1.36 

19 

0.1 

2169 

4.56e-4 

1.23e-3 

1 

67 

0.08 

2.55 

32 

0.2 

559 

3.03e-3 

3.68e-3 

51 

62 

0.46 

0.70 

1.5 

0.2 

1129 

1.33e-4 

1.98e-3 

15 

62 

0.25 

1.24 

5.0 

0.2 

2169 

7.28e-4 

1.41e-3 

20 

62 

0.61 

2.54 

4.2 


Table 2: Herschel-Bulkley flow through the cylindrical pipe, a = 1.75. 


analytical solution. To verify that Algorithms TRS and ALG2 compute accu¬ 
rate solutions, we compute a relative error \y — ya\l\ya\ between the converged 
solution y and the analytical solution ija evaluated on the vertices of the mesh. 
Furthermore, we provide the number of iterations required to satisfy each stop¬ 
ping criterion, the CPU times needed for each algorithm to terminate and the 
speedup as the ratio of the times for ALG2 and TRS. 

Looking at the relative errors in Tables [T] to O we conclude that both al¬ 
gorithms compute accurate approximations to the analytical solution. Errors 
tend to decrease as the mesh is refined. This behaviour indicates convergence 
of the numerical solutions to the exact solution as /i 0. Our results are also 
consistent with the data in [2T| . 


To 

UN 

Relative error 
TRS ALG2 

Iterations 
TRS ALG2 

CPU time (s) 
TRS ALG2 

Speedup 

0.1 

559 

3.38e-3 

3.84e-3 

20 

54 

0.20 

1.82 

9.1 

0.1 

1129 

1.50e-3 

2.01e-3 

3 

54 

0.08 

9.87 

120 

0.1 

2169 

7.87e-4 

1.32e-3 

3 

54 

0.12 

23.02 

190 

0.2 

559 

5.68e-3 

6.05e-3 

37 

43 

0.34 

4.88 

14 

0.2 

1129 

2.69e-3 

3.06e-3 

11 

43 

0.17 

13.41 

79 

0.2 

2169 

1.46e-3 

1.81e-3 

12 

43 

0.31 

27.89 

90 


Table 3: Herschel-Bulkley flow through the cylindrical pipe, a = 1.5. 
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Moreover, it takes consistently fewer iterations for Algorithm TRS to con¬ 
verge compared to Algorithm ALG2. With only one single exception for Bing¬ 
ham flow, the runtime for Algorithm TRS amounts to just a small fraction of 
the time Algorithm ALG2 takes to converge. This enormous difference in per¬ 
formance of the two methods increases even further in favour of the trust-region 
SQP algorithm when a is decreased. Overall, Algorithm TRS is hardly affected 
by changes in the power-law exponent a, neither in terms of iterations, nor 
in terms of GPU times. In sharp contrast, the number of iterations decreases 
slightly for Algorithm ALG2 as a decreases, while the computational expen¬ 
diture per iteration increases significantly. This is due to the fact that in the 
Bingham case, the auxiliary variable g is obtained through a simple function 
evaluation. With general Herschel-Bulkley fluids, however, Newton’s method is 
required to solve a nonlinear equation for g in the region of the domain where 
the fluid has yielded. The more a deviates from the value 2, the more iterations 
of Newton’s method are needed in every single step of Algorithm ALG2 to find 
the next iterate of g. 

An interesting feature of the augmented Lagrangian method is its remarkable 
mesh independence. For all combinations of the yield stress tq and the exponent 
a chosen here, the number of iterations is exactly the same for all three grids. 
The trust-region SQP method exhibits a trend of a decreasing number of iter¬ 
ations as the mesh is refined, and in some cases even terminates after the first 
iteration. With a simple calculation, one verifies for the cylindrical duct that 
one solution of the continuous problem (13.51) is given by t := —V(—A)^^/. Al¬ 
gorithm TRS is initialised with r = A~^D ^which corresponds 
to the finite element discretisation of this solution. Hence, with h sufficiently 
small, this initial projection already solves the discrete optimality conditions 
(|4.4I) and Algorithm TRS terminates. 

Overall we conclude that both methods are suited to compute accurate ap¬ 
proximations to the cylindrical duct flow. However, Algorithm TRS generally 
converges several times faster than Algorithm ALG2. 

6.2 Flow through a kiwi shaped duct 

One may argue that the advantageous performance of Algorithm TRS could, 
at least partially, be an artefact of the symmetry of the cylindrical duct. It is 
this symmetry that allows the problem to be reduced to one dimension and to 
solve it analytically. In order to eliminate such effects, we now consider a more 
complicated geometry with a non-symmetric domain. 

As shown in Figure [21 we now investigate Herschel-Bulkley flow through a 
pipe with the shape of a kiwi. Figure [3] depicts a typical solution of viscoplastic 
flow through a pipe with this cross-section and no-slip boundary conditions. In 
the thin regions of the domain around the beak and the foot, the fluid remains 
stuck and does not flow at all. Similarly, another black area in the centre of the 
body indicates how a column of fluid moves rigidly through the duct. 

For the particular choice of parameters in Figures [3] and HI both velocity 
profiles appear to be identical by visual inspection. There is a slight discrepancy 
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0.45 



0.25 


0.15 


0.05 


Figure 1: Typical velocity profile for viscoplastic duct flow, computed with 
Algorithm TRS. Here, a = 1.75 and tq = 0.2. On triangles Ti that are coloured 
black, \Ti\ < Tq. Otherwise, the colours represent \fi\. 



Figure 2: Geometry of the kiwi shaped pipe with a sample mesh (tin = 2101). 
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0.45 



Figure 3: TRS approximation of Bingham flow through the kiwi shaped pipe, 
a = 2, To = 0.2, tt-n = 2101. As in Figure [U the colours indicate locations of 
stagnant flow and represent the stress magnitude in yielded regions. 



Figure 4: ALG2 approximation of Bingham flow through the kiwi shaped pipe, 
a = 2, To = 0.2, riN = 2101. 
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T-O 

UN 

Relative difference 

Iterations 
TRS ALG2 

CPU time (s) 
TRS ALG2 

Speedup 

0.1 

1125 

9.67e-4 

14 

73 

0.20 

0.53 

2.7 

0.1 

2101 

9.67e-4 

13 

73 

0.29 

0.96 

3.3 

0.1 

4308 

9.69e-4 

8 

73 

0.34 

2.16 

6.4 

0.2 

1125 

1.08e-3 

80 

73 

1.00 

0.55 

0.55 

0.2 

2101 

1.33e-3 

27 

73 

0.50 

0.94 

1.9 

0.2 

4308 

1.19e-3 

37 

73 

1.29 

2.15 

1.7 


Table 4: Bingham flow through the kiwi pipe, a = 2. 


in the approximation of the solid regions: in the graph for Algorithm TRS, 
single triangles are coloured black unlike in the approximate solution generated 
by Algorithm ALG2. This insignificant difference is also illustrated in Figure 
m The surface in the latter graph represents the absolute error Ij/trs — yALG2|, 
the difference in the solutions obtained by the two methods. As expected, the 
difference reaches local highs near interfaces between viscous and plastic regions. 
In those plastic regions, the functional on which the augmented Lagrangian 
method is based is nonsmotth. Near the same interfaces, the gradients and 
Hessians which occur in the trust-region SQP algorithm become singular. As a 
result, the approximate solutions differ most in their prediction of the plug flow 
velocity. 

In Tables [4] to El we present similar data for the simulations like before in 
Tables [I] to El Since an analytical solution is unavailable, we compute a relative 
difference Ij/trs — yALG2|/|yALG2| instead of the relative error. 

Even though in none of the examples Algorithm TRS terminates after the 
first iteration, its clearly superior performance persists even in this more complex 
geometry. The augmented Lagrangian method still exhibits mesh independence 
to a great extent. The number of iterations for the trust-region SQP method 
tends to decrease when the mesh is refined or when tq is decreased. No obvious 
correlation can be found between iterations of Algorithm TRS and the exponent 
a. 

Again, there is only one exception to the overall trend that the computa¬ 
tion with Algorithm TRS terminates significantly faster than Algorithm ALG2. 
With a = 1.75 and even more with a = 1.5, the computational cost of the aug¬ 
mented Lagrangian method becomes excessively high. Gomparing the speedup 
factors for the two geometries, it appears that the symmetry of the problem 
is exploited more efficiently in the trust-region SQP method. However, even 
without any symmetry in the domain, we still achieve major gains in the com¬ 
putational performance by using our new method TRS. 


7 Outlook and Conclusions 

So far, we have successfully developed a new approach to model and simulate 
viscoplastic flow. Analytically, the proposed dual formulation is fully equivalent 
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Figure 5: Absolute error between the approximations of Figures |3] and IH Black 
triangles are identified as unyielded by Algorithm TRS, but not by Algorithm 
ALG2. White triangles are classified equally. 


'To 

riN 

Relative difference 

Iterations 
TRS ALG2 

GPU time (s) 
TRS ALG2 

Speedup 

0.1 

1125 

7.28e-4 

9 

64 

0.16 

1.10 

6.9 

0.1 

2101 

7.38e-4 

20 

64 

0.45 

1.81 

4.0 

0.1 

4308 

7.42e-4 

14 

64 

0.60 

4.01 

6.7 

0.2 

1125 

9.09e-4 

33 

58 

0.49 

0.92 

1.9 

0.2 

2101 

7.92e-4 

42 

58 

0.75 

0.15 

2.0 

0.2 

4308 

1.46e-3 

19 

58 

0.68 

3.34 

4.9 


Table 5: Herschel-Bulkley flow through the kiwi pipe, a = 1.75. 


To 

tin 

Relative difference 

Iterations 
TRS ALG2 

GPU time (s) 
TRS ALG2 

Speedup 

0.1 

1125 

3.63e-4 

31 

47 

0.46 

11.95 

26 

0.1 

2101 

4.13e-4 

35 

47 

0.67 

27.47 

41 

0.1 

4308 

6.47e-4 

17 

54 

0.66 

92.89 

140 

0.2 

1125 

3.21e-4 

62 

72 

0.79 

13.55 

17 

0.2 

2101 

5.46e-4 

48 

56 

0.84 

22.95 

27 

0.2 

4308 

l.lOe-3 

27 

64 

1.02 

116.62 

114 


Table 6: Herschel-Bulkley flow through the kiwi pipe, a = 1.5. 
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to the conventional definition in terms of a variational inequality, or nonsmooth 
minimisation problem. 

From a numerical perspective, however, the feature of higher regularity of the 
problem is very desirable. With Algorithm TRS, we have designed a powerful 
alternative to the well-known Algorithm ALG2. The latter method suffers from 
impaired performance due to the lack of continuous derivatives, the inability 
to exploit extra regularity in the case of shear-thinning fluids and the need to 
solve potentially large nonlinear systems repeatedly. Algorithm TRS, however, 
displays superior performance in particular for shear-thinning Herschel-Bulkley 
fluids. For our numerical examples, we generally observe large speedup factors, 
depending on the flow domain. 

The extension of Algorithm TRS to time-dependent or fully two-dimensional 
flow is straightforward. While this will affect the explicit operators and function 
spaces involved, the overall structure of the problem remains preserved. 

As with the augmented Lagrangian approach, the governing equations of 
viscoplastic flow no longer correspond to the optimality conditions of a minimi¬ 
sation problem when certain additional terms are considered in the model. This 
includes for instance finite slip (Robin) boundary conditions, stick-slip boundary 
conditions m or convection. The key question for a generalisation of Algorithm 
TRS to include those features will be, how the trust radius can be updated when 
no objective is available. In analogy to our successful application of Algorithm 
TRS to flows of nonlinear shear-thinning fluids, the trust-region SQP method re¬ 
quires no costly inner loop to handle further nonlinearites. Hence, also for more 
complex flow problems, we expect significantly increased performance compared 
to the augmented Lagrangian method. 
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